function [K] = prod2term_cov_fn(types,flag)
% hp(3): log input scale 1
% hp(4): log output scale 1
% hp(5): log input scale 2
% hp(6): log output scale 2


type1=types{1};
type2=types{2};

if nargin<2
    flag='plain';
end



switch flag
%     case 'deriv inputs'
%     DK1=@(as,bs) matrify(@(at,bt)...
%             gcov(type1,{T1,exp(hp(4))},at,bt),as,bs);
%     DK2=@(as,bs) matrify(@(at,bt)...
%             gcov(type2,{exp(hp(5)),exp(hp(6))},at,bt),as,bs);
%         
%     out2=@(as,bs) cellfun(@plus,DK1(as,bs),DK2(as,bs));
%         
%     DDK1=@(as,bs) matrify(@(at,bt)...
%             Hcov(type1,{T1,exp(hp(4))},at,bt),as,bs);
%     DDK2=@(as,bs) matrify(@(at,bt)...
%             Hcov(type2,{exp(hp(5)),exp(hp(6))},at,bt),as,bs);
%  
%     out3=@(as,bs) cellfun(@plus,DDK1(as,bs),DDK2(as,bs));
    case 'plain'
        K=@(hp,as,bs) matrify(@(at,bt)...
    fcov(type1,{exp(hp(3)),exp(hp(4))},at,bt).*fcov(type2,{exp(hp(5)),exp(hp(6))},at,bt),...
    as,bs);
    
    case 'deriv hyperparams'
      
    K=@(hp,as,bs) [{zeros(size(as,1),size(bs,1));...          % mean
                    zeros(size(as,1),size(bs,1))};...         % logNoiseSD 
                    cellfun(@(x) x.*matrify(@(at,bt)...
                        fcov(type2,{exp(hp(5)),exp(hp(6))},at,bt),as,bs),[...
                    matrify(@(at,bt)...
                        gTcov(type1,{exp(hp(3)),exp(hp(4))},at,bt),as,bs);...
                    matrify(@(at,bt)...
                        2*fcov(type1,{exp(hp(3)),exp(hp(4))},at,bt),as,bs)],...
                        'UniformOutput',false);
                    cellfun(@(x) x.*matrify(@(at,bt)...
                        fcov(type1,{exp(hp(3)),exp(hp(4))},at,bt),as,bs),[...
                    matrify(@(at,bt)...
                        gTcov(type2,{exp(hp(5)),exp(hp(6))},at,bt),as,bs);...
                    matrify(@(at,bt)...
                        2*fcov(type1,{exp(hp(3)),exp(hp(4))},at,bt),as,bs)],...
                        'UniformOutput',false)];
    case 'vector'
        K = @(hp,as,bs) fcov(type1,{exp(hp(3)),exp(hp(4))},as,bs,'vector')...
                .*fcov(type2,{exp(hp(5)),exp(hp(6))},as,bs,'vector');
                    
end
